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CMBEASY:: an Object Oriented Code for the Cosmic Microwave Background 
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We have ported the CMBFAST package to the CH — h programming language to produce CMBEASY, 
an object oriented code for the cosmic microwave background. The code is available at 
www.cmbeasy.org. We sketch the design of the new code, emphasizing the benefits of object orienta- 
tion in cosmology, which allow for simple substitution of different cosmological models and gauges. 
Both gauge invariant perturbations and quintessence support has been added to the code. For ease 
of use, as well as for instruction, a graphical user interface is available. 



I. INTRODUCTION 

The CMBFAST computer program for calculating cos- 
mic microwave background (CMB) temperature and po- 
larization anisotropy spectra, implementing the fast line- 
of-sight integration method, has been publicly available 
since 1996 It has been widely used to calculate spec- 
tra for open, closed, and flat cosmological models con- 
taining baryons and photons, cold dark matter, massless 
and massive neutrinos, and a cosmological constant. It 
is a very well-tested program that has enabled many cos- 
mologists to compare models of the Universe against the 
experimental CMB data. 

However, from the point of view of code design, there 
are few programs that could not be improved. This is 
also true for CMBFAST: it is a rather monolithic code 
that is quite difficult to oversee and modify. For exam- 
ple, variations in the cosmic expansion history due to 
dark energy necessitate similar changes to the code in 
numerous separate locations. 

In order to address these shortcomings and simplify 
modifications of the code - in our case the implementa- 
tion of quintessence models and gauge invariant variables 
- we have ported the CMBFAST package to the C-I--I- pro- 
gramming language. The C-I--I- language is object ori- 
ented and it turns out that to think in objects (more of 
this soon), is very advantageous in cosmology. 

The program has not been rewritten from scratch, but 
redesigned step by step. Some people may argue that it 
is hence not independent, i.e. some unknown errors and 
limitations in CMBFAST could be present in the new code. 
The object oriented modular design, however, ensures 
that each part of the code is independently testable. If, 
for instance, one does not trust the integrator, one can 
use another one to check it, without changing anything 
else in the package. In addition most of the lines in the 
code have been rewritten, to benefit from the redesign. 

CMBEASY calculates the temperature and polarization 
anisotropies for spatially flat universes, containing dark 
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energy, dark matter, baryons and neutrinos. The dark 
energy can be a cosmological constant or quintessence. 
There are three scalar field field models implemented, or 
alternatively one may specify an equation-of-state his- 
tory. Furthermore, the modular code is easily adapted 
to new problems in cosmology. 

In the next section lll Al we explain object oriented pro- 
gramming and its use in cosmology. The design of CM- 
BEASY is presented in section IIIBI while the graphical 
user interface is discussed in section III CI Computing 
precision is discussed in section lll Dl while the documen- 
tation for CMBEASY is introduced in III El For an intro- 
duction to perturbation theory and CMB physics using 
the conventions of CMBEASY, as well as a summary of 
the evolution equations, please see the documentation 
on www.cmbeasy.org. 



II. THE CODE 

We briefly introduce the concept of object oriented pro- 
gramming before turning to the design of CMBEASY. We 
do this using examples from the code. 



A. Objects 

In C++, a class is a user-defined data type As 
with any data type, there can be many variables with 
the same type, e.g. floats, arrays, complex numbers. In 
the case of user-deflned types, the speciflc variables are 
called objects. 



1. Encapsulation 

Quite often, some data and functions acting on the 
data are so tightly connected, that it is sensible to think 
of them as one object. As an example, let us discuss 
splines. Given a discrete set of n points Xi with Xi < 
Xi-fi and corresponding f{xi) = yi, a spline can smoothly 
interpolate, i.e. give f{x) for any x G [a;o,a;„]. Provided 
that the sampling is dense enough, arbitrary functions 
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FIG. 1: The visibility g = «;exp(K(r) — fi^(T"o)) as a function 
of conformal time r in Mpc. Its peak at about r ~ 300 Mpc 
defines the epoch of last scattering. Before the visibility func- 
tion peaks, photons are very likely to scatter again until the 
Universe becomes translucent. After the peak, photons do 
not scatter at a substantial rate. It is thus the balance be- 
tween frequent scattering and sufficiently low optical depth 
that will give the largest contribution towards the anisotropy 
today. And this fact is exactly encoded in g. 



may be described by a spline for all practical purposes. 
This is widely used in CMBFAST. For instance, the Cj's 
are calculated only every 50 1-values for / > 200. As the 
spectrum is very smooth, this still gives a precise result. 

Now, the function describing the photon visibility (see 
also Figure^, calculated in the part of the code that 
tracks the thermal history of the Universe, can be used 
to define a spline. Without object orientation, one would 
need to keep track of various variables, most notably ar- 
rays for the X, y data and derivatives needed for spline in- 
terpolation. Also, in order to assure quick access within 
the spline data table, one either needs to know the precise 
layout of the data arrays (cmbfast does this), or even 
more variables (storing for instance the last interpolation 
X value) would be necessary. In total, this sums up to a 
lot of bookkeeping for a conceptually simple entity like a 
spline. 

Alternatively, one may define a class holding all the 
necessary variables a spline needs, together with defi- 
nitions of an interface with which other parts of the 
program can access and manipulate the spline data. 
An object behaves as described by the corresponding 
class. There can be an arbitrary number of objects 
of a certain class (just like there is one floating point 
type float, but many variables of type float in a pro- 



gram). -"^ The class (in our case) called Spline, can hence 
be viewed as yet another data type, with no more book- 
keeping needed than say for a floating point number. 
To illustrate this, let us discuss the visibility function 
g = Kexp(K(r) — ^(ro)), where k[t) is the optical depth. 
Its typical shape is depicted in Figure^ and its peak de- 
fines the epoch of last scattering. As soon as the Spline 
called visibility has been given the data, its maximum 
can be determined by a single line of code: 

tau_ls = visibility .maximumO ; 

Here, the function maximum() acts on the visibility. The 
important point to notice is that all functions defined in 
the Spline class arc immediately available to everyone 
who uses Splines. So, whenever one needs to find the 
maximum, integrate a spline, calculate the convolution 
of two Splines etc, this can be done using very few lines 
of code: the functionality is fully encapsulated in the im- 
plementation of the Spline class. Any increase in perfor- 
mance or sophistication of the Spline class immediately 
translates over to all Splines used in the program. 



2. Inheritance 

Tightly connected to the fact that data and methods 
are combined within one object, is the concept of in- 
heritance, which proves very powerful in cosmology. A 
class can inherit from another class (in this context called 
base class). All variables and the full functionality that 
the base class implements arc instantly available to the 
inheriting class, ^ called sub-class. The sub-class can then 
re-implement functions of the base class to provide a dif- 
ferent functionality, or add new functions and variables. 
The important point to note is that all classes deriving 
from the same base class necessarily need to provide all 
the functions the base class provides. Hence, whenever 
the code uses the base class, one can substitute an inher- 
iting class for the base class by simply changing one line 
of code. 

As an illustration, let us look at the Perturbation 
class of CMBEASY which is designed to evolve the pertur- 
bation equations for one fc-mode through conformal time 
and calculate the anisotropy sources. The Perturbation 
class defines functions to do these calculations which 



^ Wo usually denote here (and in the code) classes with capital 
first letter. In some cases where there is only one object of a 
class used in the code, we denote the object with the same name 
as the class, but with lower case initial letter. Hence, the line 

Cosmos cosmos; 

creates an object 'cosmos' of the class 'Cosmos'. 
^ This is as if a child was born with the whole knowledge of its 
parents. No training and learning would be necessary. It could 
instantly go and increase its capabilities starting from the level 
of its parents. 
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other parts of the program can be sure to find imple- 
mented in all sub-classes. In practice, there are four 
classes that inherit from it, for perturbations in gauge- 
invariant variables and in synchronous gauge both with 
and without quintessence (see also Figure 01 ■ From the 
point of view of the rest of the program, all of them are 
equally well suited.^ Hence, when given a Perturbation 
object, one may e.g. ask it to initialize the perturbation 
variables: 

perturbation->initialScalarPerturbations ; 

Yet, depending on the specific sub-class the object be- 
longs to ^, the implementation of initialScalarPerturba- 
tions() and hence the code executed will be different. 
Technically this is called polymorphism. 



B. Design 

There are roughly three main steps needed to cal- 
culate the CMB anisotropy spectrum within the line 
of sight strategy: (1) solve the expansion and thermal 
background evolution, (2) propagate the perturbation 
equations in Fourier space, (3) map the temperature 
anisotropy onto the sky today. In each step, several 
classes work together. To illustrate the hierarchy of the 
most important classes, an overview is given in Figure |21 
The core part of the package is the CmbCalc class. It is 
the scheduler calling other objects to complete all three 
steps. So let us turn to step (1), the expansion history 
of the universe. 



1. The background evolution 

The Cosmos class is the central instance providing 
background quantities like the energy density p{t) of all 
species etc. This centralization of the background evolu- 
tion within the Cosmos class facilitates the modification 
of the code significantly. A different background cosmol- 
ogy (such as quintessence) can be implemented by just 
inheriting from Cosmos and re-implementing the expan- 
sion history part of the code.^ For most other objects in 
the code, it makes no difference exactly what caused the 
universe to expand, as long as they can rely on Cosmos 
or any sub-class to provide access to the expansion his- 
tory. Both the expansion and the thermal history are 
computed at the beginning of the calculation and stored 
in splines. 



^ Except for the fact that if one wants quintessence, the perturba- 
tion class should of course support it. 
* Representing the gauge choice and quintessence support 
^ All in all 800 lines of a total of 2500 lines of Cosmos. 
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2.. Propagation in k-space 

All fluctuation equations are handled by 
Perturbation and its sub-classes. Its most promi- 
nent function is propagateScalar (ti, T2) which 
evolves the perturbations from conformal time ti to T2. 
For numerical stability, one usually treats the era of 
tight coupling between baryons and photons differently. 
However, this results in a non-analytic jump in the r.h.s 
of the perturbation equations at Ttc- Some integrators 
for ordinary differential equations get confused by this. 
Therefore, propagateScalar () will determine if the 
jump occurs within [ti, and if so, integrate up to Ttc, 
make a tiny^ jump and integrate on till T2. At the end 
of each propagation, the sources of cmb anisotropics are 
calculated by the scalarSources () frmction. 

3. The Cl's 

The convolution of the sources with the Bessel func- 
tions, as well as the final fc-integration are performed 
in the Scalarlntegrator sub-class. In contrast to CMB- 
FAST, we formulate the convolution with the j'^s as an or- 
dinary differential equation (ODE). In fact, the Spline 
class provides the convolution functionality. If higher 
precision in this step is wanted, one may simply in- 
crease the desired precision for Spline, which will pass 
this through to the integrator of ODE's. As far as the 
fc-integration is concerned, we stick to the interpolation 
points of CMBFAST. Wc have written two other imple- 
mentations of Scalarlntegrator, although none could 
compete with CMBFASt's speed. Furthermore, no gain in 
accuracy justifying the speed tradeoff was found. As far 
as the fc-integration is concerned, a fast and independent 
realization would be desirable. With the high level fea- 
tures of C++, algorithms that seem hard to implement 
in other languages may be feasible. 



4- Little helpers 

The most valuable class assisting the computation may 
be Spline. As splines are used all over the code, and 
quite often in time-critical environments, the main design 
goal has been speed. In order to perform a spline inter- 
polation at X, one needs to find the interval 
While CMBFAST hard-wires the size and layout of the Xi- 
data, this has not been an option for a versatile spline 
implementation. Yet, in all time critical situations within 
our calculation, two facts come to the rescue: Firstly, in- 
terpolations occur in order and close to the last interpo- 
lation. Hence, the next value of x is likely to lie within 



the same interval (or close to it). Secondly, many splines 
with exactly the same x-data are evaluated at the same 
X. For instance, the r.h.s of the perturbation equations 
need background quantities all at the same r. As all 
background quantities get their data at the same Ti in 
Cosmos, they necessarily share the same x-data. There- 
fore, some parts of the interpolation formulae have to be 
computed only once^ for each t. In addition, finding the 
position of the interpolation interval is simple: it will be 
the same for the same r. 

To accommodate for these groups of splines and 
caching, splines come in two flavors: with and without its 
own x-data. While a spline with its own x-data can live 
on its own, a spline without x-data can only be created if 
a "mother" spline providing the x-data is given upon the 
creation of the "child" spline. Within these "families" of 
shared x-data, interpolation at the same x is particularly 
fast. 

Not only must a spline be fast, it must also be robust. 
Hence, a spline will allocate more memory if needed®, 
check that no access to data occurs outside its bound, 
check that all data is ordered Xi < Xi^i and make sure 
that the interpolation tables have been built before any 
interpolation takes place. Splines can easily be visual- 
ized. To write the data^ to a file, simply call 

myspline . dump ( " f ilename " ) 

The anisotropy sources, are evaluated on a {fc, r} grid. 
For the final r and fc integration, one needs to interpolate 
at intermediate points. This is conveniently done using 
the class SplineWeb. It is essentially an array of Splines, 
one for each fc and one for each t value. If one needs to 
interpolate at a given k value, one can ask SplineWeb to 
return a Spline with x-data in the t direction: 

spline_at_k = web . createAlongY(k) ; 

This is used to get the splines which are convoluted with 
the Bessel functions. It is also very convenient for mak- 
ing slices through the cold dark matter power spectrum 
(at either fixed fc or r). Given for instance the cold dark 
matter SplineWeb power, one can dump it to a file: 

power .dumpAlongXCr, "filename") ; 

The above command would save the spectrum at arbi- 
trary conformal time r to disk. 

There are some more of these "helpers" , most of them 
small classes that sometimes make life much easier. The 
last one we will introduce here is Anchor which keeps 
track of objects that are dynamically created (as are 
many splines). The objects it is asked to keep track of 



' With "once" , we mean until interpolation is requested for a dif- 
ferent T. 

More specific, it adds Ar, such that for double precision variables * and free unused memory 

Ttc + = Ttc- ® and interpolated values, if desired 
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FIG. 3: Graphical user interface (GUI) for CMBEASY. 



are deleted by it on demand, or when the Anchor object 
itself is deleted. Therefore, it greatly simplifies memory 
handling and prevents memory leaks. 



5. Quintessence Implementation 

The different background evolution of quintessence 0, 
scenarios is implemented using the QuintCosmos and 
the Quintessence class. Each subclass of Quintessence 
corresponds to a certain model, such as the exponen- 
tial potential 0, 0, 0], inverse power law leaping 
kinetic term Q etc. Certainly, a more monolithic de- 
sign with the quintessence models implemented in the 
QuintCosmos class would have been possible. However, 
we believe that the details of a model are best kept to a 
class of its own. For instance tuning model parameters 
in order to get the right amount of 17q is different for 
each model and a monolithic design would have to call 
differently named functions for different models. Using 
sub-classing, QuintCosmos (and Perturbation) always 
calls functions with the same argument and name for 
all models. Yet, as the object implementing the func- 
tion differs for each model, the code executed by call- 
ing the function can be totally different. Thus, a new 
quintessence model can be implemented by simply sub- 
classing Quintessence and implementing a minimal set 
of functions, such as one for the scalar field potential. 



C. Graphical User Interface 



ily, there is the very sophisticated and publicly available 
'Qt' library @ with which the creation of a GUI is fa- 
cilitated. Its object oriented C+-I- design is a perfect 
match for the CMBEASY package. There is therefore an 
executable program called 'cmbeasy' giving interactive 
access to almost the full capabilities of the package, in- 
cluding quintessence. As is seen from Figurc|31 the spec- 
tra arc visualized in separate plots arranged in a so called 
Tab- Widget. One can for instance zoom in, select and 
save curves or print the plot. 



D. Computing precision 

Comparing the results from CMBFAST to those ob- 
tained from CMBEASY, one generally notes accordance 
to better than 0.5% (see also Figure 0J. 

In principle, the line of sight algorithm is equivalent to 
evolving a large hierarchy of multipoles. In practice, as 
one truncates the hierarchy in the line of sight approach 
at some low Z, there is a reflection of power after the 
higher multipoles gradually build up. However, as by 
far the largest contribution towards the anisotropy oc- 
curs around recombination (when the higher multipoles 
are still small), this does not spoil the result dramati- 
cally. Yet, to check this, one may increase the number of 
multipoles involved. This is true both for CMBFAST and 

CMBEASY. 

Within CMBEASY, there are five instances in which the 
precision may be increased. In three of the five cases, 
changing one variable suffices. The first case is the back- 
ground evolution. Here one may want to increase (or 
decrease) the precision requested from the ODE integra- 
tor. During perturbation evolution in fc-space, the preci- 
sion requirement given to Perturbation can be changed 
in a similar manner. The same is true for the convo- 
lution of the sources with the Bessel hmctions, taking 
place in Scalarlntegrator. Finally, and a bit more 
difficult, one may change the number and position of 
fc-values contributing towards the fc-integration. This 
can be done by modifying Integrator, the base class 
of Scalarlntegrator. Alternatively, one could devise a 
different implementation for the /c-integration^^. Finally, 
one may increase the number of Z's at which the C/'s are 
calculated. 



E. Documentation 

Using the doxygen program, the documentation is 
automatically generated from the source code of the CM- 
BEASY package. In its html version, it is interactively 



For educational purposes and also to simplify the pa- 
rameter input and subsequent visualization of results, a 
graphical user interface (GUI) is of great value. Luck- 



A widget is a part of a user interface that can interact. Examples 
are buttons, sliders, etc. 
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a new sub-class of Integrator 



6 



0.3 
0.2 
g 0.1 


U- -0.1 
-0.2 



- 






1 

1 \ 

1 \ f\ ' 
















Ma/ V 

jr"^ ..1 














1 



80 
70 
60 
50 
40 
30 
20 



(N 



u 
+ 



500 



1000 



1500 



FIG. 4: Temperature anisotropy spectrum for h = 0.65, Qq — 
0.6, floh' = 0.02, = 1 - - njj obtained from CMBFAST. 
The relative deviation ACi/Ci of CMBEASy's synchronous 
(long dashed line) and gauge invariant (solid line) solution 
with respect to the original CMBFAST spectrum are also given. 
The accordance of all spectra is always better than 1%. In 
the gauge invariant case, both the background and perturba- 
tion evolution as well as the Ci integration are entirely inde- 
pendent of the CMBFAST code. However, they use the same 
thermal history algorithm that should in principle be inde- 
pendently implemented for cross checks. 



navigable and includes the full source code. Due to the 
automatic generation, the documentation and the code 
are naturally synchronized. Instructions for setting up 
and running the program are available in this html doc- 
umentation. In addition, a postscript version of the doc- 
umentation is generated. 



III. CONCLUSIONS 

We have devised an object oriented code for calculating 
cosmic microwave background anisotropics. The main 
benefits of this code derive from its modularity. 

The modular design allows one to test and refine each 
part of CMBEASY individually. Different modules, for 
different numerical methods and theoretical approaches, 
can easily be switched. This feature should prove ex- 
tremely useful in the development and refinement of the 



code. As cosmological theory reaches for a new level of 
precision to match the new generation of experiments, 
MAP 0, PLANCK ^3 and other experiments will soon 
put our predictions and hypotheses to the test, so future 
work will seek to improve the accuracy of the code. (At 
the present time, CMBEASY and the standard CMBFAST 
arc generally in agreement to better than 0.5%.) 

The modular design allows one to more easily develop 
and build upon CMBEASY, to implement new cosmolog- 
ical models and include more physics. As our theoret- 
ical cosmological models are a work-in-progress — as 
CMB data comes in — so are our tools for modeling 
the Universe. And so it makes sense to use a program 
design which makes this development easy. Hence, the 
bmctional blocks used in CMBEASY considerably reduce 
the time to implement new cosmological models. For 
most purposes, a new model can be achieved by sub- 
classing from existing solutions, thus eliminating the need 
to rewrite large parts of the code. Furthermore, new 
modules can be easily introduced to carry out calcula- 
tions that rely on the linear perturbation spectrum. Ex- 
amples include lensing, secondary sources in the CMB, 
reionization, and cross-correlation with other tracers of 
large scale structure. 

An additional benefit of the CMBEASY package is the 
graphical front-end. As we have described, the GUI al- 
lows a user to visually dial through a range of parameters. 
We expect that this functionality will be useful for devel- 
oping understanding and intuition for CMB physics and 
cosmological models. 

As stated at the outset, there is, perhaps, no program 
design that could not be improved. However, we hope 
that CMBEASY will serve the cosmology community well. 
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